#PBS -N 04-01AB_star_qualimap
#PBS -A GT-wratcliff3
#PBS -q inferno
#PBS -j oe
#PBS -o logs/04-01AB_star_qualimap.out
    
# change to project directory
PROJECT_DIR=$PBS_O_WORKDIR
echo "Changing directory to ${PROJECT_DIR}…"
cd $PROJECT_DIR
     
# activate environment
PIPELINE_ENV=`basename scripts/*.yml .yml`
echo "Activating conda environment ${PIPELINE_ENV}…"
source activate $PIPELINE_ENV
    
# extract sample name
samplenames=(`cat raw_data/samplenames.txt`)
sampleindex=$PBS_ARRAYID
samplename=${samplenames[$(($sampleindex-1))]}
     
echo -e "\n########## Start processing sample ${sampleindex}: ${samplename} ##########\n"

echo -e "\n########## Starting STAR ##########\n"

# ref
star_genome_index=/storage/home/hcoda1/6/ktong34/p-wratcliff3-0/ngs_ref/ref_index/ensembl_release-103_Scer/star_v2.7.8a_genome_index/read-length-75

# input/output
in=results/03A_sortmerna/out
in_R1=$in/${samplename}_rRNA_filtered_R1.fastq.gz
in_R2=$in/${samplename}_rRNA_filtered_R2.fastq.gz
out=results/04-01A_star
out_prefix=$out/${samplename}_
mkdir -p $out

# parameters
intron_min=4   # yeast
intron_max=5000   # yeast

# main
STAR \
--readFilesIn $in_R1 $in_R2 \
--readFilesCommand zcat \
--genomeDir $star_genome_index \
--outFileNamePrefix $out_prefix \
--runThreadN $PBS_NP \
--outSAMtype BAM SortedByCoordinate \
--outBAMsortingThreadN $PBS_NP \
--outSAMunmapped Within \
--outSAMattributes Standard \
--limitBAMsortRAM 3000000000 \
--outBAMsortingBinsN 150 \
--alignIntronMin $intron_min --alignIntronMax $intron_max

echo -e "\n########## Starting Qualimap ##########\n"

# ref
gtf=/storage/home/hcoda1/6/ktong34/p-wratcliff3-0/ngs_ref/ref_data/ensembl_release-103_Scer/gtf/Saccharomyces_cerevisiae.R64-1-1.103.gtf

# input/output
in=results/04-01A_star/${samplename}_Aligned.sortedByCoord.out.bam
out=results/04-01B_qualimap/$samplename  # save to one directory for each run otherwise overwrite 
mkdir -p $out

# main
unset DISPLAY  # do not use Qualimap GUI
qualimap rnaseq \
-bam $in \
-gtf $gtf \
-outdir $out \
-p strand-specific-reverse \
-a proportional \
-pe \
--java-mem-size=8G

